1 Estacionariedade

1.1 Introdução

A série temporal tipicamente faz a relação de uma entidade em vários períodos de tempo. Por exemplo, temos a serie temporal AirPassengers (embutida no proprio R) que nos mostra o total mensal de passageiros de voos comerciais dos EUA durante 1946-1960

data(AirPassengers)
AP <- AirPassengers
plot(AP)

Em geral, as primeiras análises procuram observar padrões no tempo, como a tendência de crescimento ou decrescimento, ciclos de altas(baixas), ou sazonalidade (padrões repetitivos no curto prazo) e eventuais irregularidades. Seja uma relação de séries temporais $y$ e $x_k$:

\[{y_t} = {b_0} + {b_1}{x_{1t}} + \ldots + {b_k}{x_{kt}} + {u_t}\]

As séries temporais possuem uma ordenação explícita no tempo. As relações podem ser contemporâneas (mesmo tempo t) ou dinâmicas (tempos distintos):

${y_t} = {b_0} + {b_1}{x_{1t}} + \ldots + {b_k}{x_{kt}} + {u_t}$ contemporâneas

${y_t} = {a_0} + {d_0}{z_t} + {d_1}{z_{t - 1}} + {d_2}{z_{t - 2}} + {u_t}$ dinâmicas ou defasadas

O exemplo anterior é típico de defasagem distribuída de ordem 2. O processo genérico poderia ser de q defasagens de z, ou ordem q , dado por:

${y_t} = {a_0} + {d_0}{z_t} + {d_1}{z_{t - 1}} + \ldots + {d_q}{z_{t - q}} + {u_t}$

1.1.1 Passeio Aleatório (Random Walk)

O passeio ou caminho aleatório é um processo $\{X_t\}$ tal que depende apenas de seus valores passados e de um processo puramente aleatório $\{ε_t\}$ com média $\mu$ e variância $σ^2_ε$ tal que:

$X_t = X_{t-1} + ε_t$ caminho aleatório sem intercepto
$X_t = \delta + X_{t-1} + ε_t$ caminho aleatório com intercepto

Substituindo a expressão de $X_{t-1}$ em $X_t$ tem-se $X_t = X_{t-2} + ε_{t-1} + ε_t$. Generalizando tem-se $X_t$ como função do valor inicial $X_0=0$.

1.1.1.1 Simulação de um passeio aleatório genérico

O passeio aleatório no R é feito numa simulação, neste caso, para uma função normal genérica.

# SIMULANDO UM PASSEIO ALEATORIO GENERICO
# simulate random walk 
#
set.seed(321)
e = rnorm(250)
y.rw = cumsum(e)
ts.plot(y.rw, lwd=2, col="blue", main="Random Walk")
abline(h=0)

É possível verificar oscilações na série ao longo do tempo. A presença de eventuais tendências (se alterando ao longo de subamostras) e oscilações que indicam volatilidade ao longo das diferentes subamostras. Tais indicações, embora graficamente ainda representem apenas indicações, darão possíveis interpretações de não-estacionariedade da série.

1.1.1.2 Exemplo para a série do PIB Brasileiro

Neste exemplo utiliza-se a série n. 4380 do PIB mensal a valores correntes de 01/1990 até 04/2023

PIB <- Quandl('BCB/4380')
head (PIB)

##         Date    Value
## 1 2023-04-30 881892.4
## 2 2023-03-31 900662.6
## 3 2023-02-28 823682.7
## 4 2023-01-31 804545.5
## 5 2022-12-31 868939.6
## 6 2022-11-30 859384.4

Aqui podemos ver que o data frame se inicia em 2023 e vai indo ate 1990, isso se tornara um problema para criarmos uma time-series no futuro, então vamos colocar ela de cabeça para baixo. Tambem vamos excluir a coluna Date, para não criar redundancia de datas na criação da TS.

PIBc <- arrange (PIB, Value) #colocando em ordem crescente (cabeça pra baixo)
PIBcpb<- PIBc [,-c(1)] #excluindo a primeira coluna (Date)
head (PIBcpb)

## [1] 0.2 0.4 0.7 0.8 0.8 0.8

AGORA SIM!

Vamos transforma-la em TS e puxar um gráfico

PIB.ts <- ts(PIBcpb, start = 1990, frequency = 12)
dygraph(PIB.ts, main = "Produto Interno Bruto (PIB) mensal - valores correntes") %>% dyRangeSelector()

Agora vamos a estacionariedade do PIB

dpib<-diff(PIB.ts,1)
ddpib<-diff(PIB.ts,2)
plot(ddpib,main="Séries de Diferenças do PIB",type = "o",col="black",lwd=2,lty=1, ylab = "Indice",xlab="")
lines(dpib,type="o",col = "red",lwd=2,lty=2)
legend("bottomleft",c("d(PIB,2)", "d(PIB)"),
       cex=0.7,lwd=2,lty=1:2,col=c(1,2))

1.1.2 Ruído Branco (White noise)

Assume-se em geral que os resíduos ut sejam bem comportados, ou seja, Ruído Branco: série de resíduos que representa um processo gaussiano, ou seja, com distribuição normal, média zero, variância constante e não-autocorrelacionados. O conceito de ruído branco será importante para compreender o conceito de estacionariedade de séries temporais.

\[ {u_t} \sim N(0,{\sigma ^2}) \]

1.1.2.1 Simulação de um processo Gaussiano White Noise

Utilizando o R/RStudio para gerar um White noise, pelo script. O resultado será como na Figura.

options(digits=4, width=70)
# simulate Gaussian White Noise process
set.seed(123)
y = rnorm(250)
ts.plot(y,main="Processo Gaussiano White Noise",
        xlab="time",ylab="y(t)",col="blue", lwd=2)
abline(h=0)

# plot equivalente usando a função plot()
plot(y, main="Processo Gaussiano White Noise", type="l",
     xlab="time",ylab="y(t)", col="blue", lwd=2)
abline(h=0)

1.1.3 Tendência determinística

Outro conceito é o de tendência (trend), que reflete as oscilações de longo prazo em uma série. No script e Figura, tem-se uma simulação de um processo de tendência determinística.

1.1.3.1 Simulação de um processo com tendência determinística

# SIMULANDO UMA TENDENCIA DETERMINÍSTICA
# 
set.seed(123)
e = rnorm(250)
y.dt = 0.2*seq(1,250) + e
ts.plot(y.dt, lwd=2, col="blue", main="Tendência Determinística + Ruído")
abline(a=0, b=0.2)

1.1.4 Estacionariedade

Em geral, a econometria tradicional (causa e efeito) pressupõe que todas as séries sejam estacionárias. Seja a seguinte série temporal $X_t$: \[ X_t = \mu + ε_t \]

Em que $ε_t$ é um ruído branco e $\mu$ é uma constante. Tem-se uma série estacionária se o valor esperado da série $X_t$ for constante no tempo, ou seja,

\[ E(X_t) = E(μ + ε_t ) = μ + E(ε_t ) = μ + 0 = μ \]

O caso da série não estacionária, por exemplo, seria o caso de uma série com tendência:

$X_t = \alpha + βt + ε_t$ em que t denota tempo;
$E(X_t) = E(\alpha + βt + ε_t) = E(\alpha + βt) + E(ε_t ) = μ + βt$ $\ne$ constante

Este é um caso de tendência estacionária, pois uma vez retirada a tendência da série, esta se tornaria estacionária. Portanto, “será uma série estacionária se os dois primeiros momentos (média e variância) forem independentes do tempo e a autocovariância for dependente apenas da janela temporal (j) entre os dados”:

$E(X_t) = μ$
$E\left[ {\left( {{X_t} - \mu } \right)\left( {{X_{t - j}} - \mu } \right)} \right] = \left\{ {\begin{array} {{\sigma^2},\;para\;\left( {j = 0} \right)}\\ {0\;\;,\;para\;\left( {j \ne 0} \right)} \end{array}} \right.$

Resumidamente, a série não-estacionária tem uma raiz unitária e representa um processo estocástico. Mas, para esclarecer, são precisos outros conceitos auxiliares antes de esclarecer o que é ter uma raiz unitária. A utilização de séries não-estacionárias poderá gerar regressões espúrias, que serão explicadas mais a frente.

1.1.5 Operadores Diferença (difference) x Defasagem (lag)

Seja a série $X_t$ em nível. O operador diferença é definido como segue:

\[ ∆X_t=X_t-X_{t-1} \]

\[ ∆^2 X_t=∆X_t-∆X_{t-1}= (X_t-X_{t-1})-(X_{t-1}-X_{t-2})= X_t-2X_{t-1}+X_{t-2} \]

O operador de defasagem (L), ou lag operator, é definido como segue:

\[ LX_t=X_{t-1} \]

\[ L^2 X_t=X_{t-2} \]

\[ L^jX_t=X_{t-j} \]

Propriedades de L:

Se c é uma constante, então: $L^j c=c$
Propriedade distributiva: $(L^i+L^j)X_t=L^iX_t+L^jX_t=X_{t-i}-X_{t-j}$
Propriedade multiplicativa: $L^i(L^j)X_t=L^{i+j}X_t=X_{t-i-j}$
Se $|a|\gt 1: (1+a^{-1}L^{-1}+a^{-2}L^{-2}+a^{-3}L^{-3}+\ldots)X_t=\frac{-aLX_t}{(1-aL)}$
Se $|a|\lt 1: (1+aL+a^2L^2+a^3L^3+\ldots)X_t=\frac{X_t}{(1-aL)}$

1.1.5.1 Operador defasagem (lag)

No R, pode-se obter a defasagem de uma série usando a função lag(x,y), em que x é a série e y é o número de períodos a defasar a série.

# fazendo lag
tail (PIB.ts)

##         Jan    Feb    Mar    Apr May Jun Jul Aug Sep Oct    Nov
## 2022                                                     853787
## 2023 859384 868940 881892 900663                               
##         Dec
## 2022 855802
## 2023

pib_3<- lag (PIB.ts, -3)
tail (pib_3)

##         Feb    Mar    Apr    May    Jun    Jul
## 2023 853787 855802 859384 868940 881892 900663

plot(tail(PIB.ts), type="o",col = "black",lwd=2,lty=1)
lines(tail(pib_3),type="o",col = "red",lwd=2,lty=2)
legend("topright",c("PIB", "PIB t-3"),lwd=2,lty=1:2,col=c(1,2))

data.lag<-cbind(PIB.ts,pib_3)
#View(data.lag)

1.1.5.2 Operador diferença (diff)

Já para fazer a diferença, faz-se, em R: diff(x, lag = 1, differences = 1).

#View(consumo.ts)
pibd<-diff(PIB.ts,1)
pibd

##          Jan     Feb     Mar     Apr     May     Jun     Jul     Aug
## 1990             0.2     0.3     0.1     0.0     0.0     0.1     0.1
## 1991     0.3     0.3     0.1     0.6     0.5     0.5     0.5     0.7
## 1992     2.6     3.2     3.3     4.0     6.5     7.4     8.6    10.6
## 1993    36.3    42.9    87.7    73.5   113.5   131.9   227.9   308.1
## 1994  1027.9  1230.2  2727.4  4308.4  7675.4 12622.4  7661.4  3048.2
## 1995  1102.4  2925.9  6396.0   375.7  1063.3   998.6  1236.2   586.5
## 1996   591.0   114.2  2444.4  1185.2  1725.8   972.9   611.2    30.2
## 1997   858.2  1372.0   382.5   528.2  1253.1  1085.2    59.6   420.4
## 1998   450.3   320.2    71.4    73.8   919.3  1681.7   474.0   182.6
## 1999   358.7     5.6   436.2  1063.6   420.4  1388.8   618.6   145.8
## 2000     3.3  3292.5  2854.6    45.5   790.6   744.5  1327.7   895.4
## 2001   673.8  1099.5   620.9   106.0   625.8   397.0  2058.1   443.3
## 2002  3074.3  1743.0   562.7  1941.2  3038.5   128.1  1585.3   573.3
## 2003  1742.4  2841.2  1959.4  2132.2    89.3   855.7   647.5   807.8
## 2004  4400.2  3595.4   759.1   855.9   146.6   560.4   971.0   924.7
## 2005   874.3  2605.9   920.6   100.1    81.9    28.4  1388.9  2156.4
## 2006  4368.0   136.8  1571.1   911.3   845.0  2316.2  2818.0  1947.2
## 2007  3591.5   118.4  1942.4  2223.7   991.5  2829.6   505.9  3381.1
## 2008  5120.3   768.3  1514.7   904.7  1554.1  5950.5  2547.5  2841.5
## 2009    99.3  2444.7   120.1   707.0   654.8  8691.4  1049.1   847.7
## 2010  4532.1  1190.3  2516.3  4296.5   599.1  2105.1  8939.2  1300.4
## 2011  2932.8  3680.1  1999.5  1809.8  5394.6   368.9  3142.5   192.1
## 2012  1714.8   389.1  5658.0  3845.3  2375.1   335.2  3768.7   813.5
## 2013  7848.4  4805.1  3224.6   109.5  6595.6  5169.7  1250.0  4842.6
## 2014  5438.5   483.9  3071.8  3932.1  1456.5  1486.4  2519.1  1553.8
## 2015    94.6   890.7   846.0  4885.8  1330.3  1242.6  1385.6   266.2
## 2016     5.2  1989.5   271.1   641.7  4301.1   567.8   280.6   283.1
## 2017  1939.5   551.2  1833.1  2291.9  1745.6  1060.7  1277.9  2570.6
## 2018  6527.3   497.4  4764.5  2039.5   368.0  1321.6   460.3  1117.4
## 2019  4606.9  2245.8   718.1  1870.9  5740.6  1409.5  3181.7   949.4
## 2020  9462.8  2370.5  5761.6   198.0    61.5  1181.0  6048.2  1815.5
## 2021 25262.4  3895.8 16314.7  1285.0   350.1  9305.8 12758.4  7245.4
## 2022  2581.8 11319.5 21693.4 15175.6  3961.6  2274.3   202.3 10687.6
## 2023  3582.7  9555.2 12952.8 18770.2                                
##          Sep     Oct     Nov     Dec
## 1990     0.1     0.3     0.3     0.1
## 1991     0.3     1.9     1.7     1.3
## 1992    14.1    18.2    21.4    17.3
## 1993   428.9   492.9   721.5   742.1
## 1994    37.2   336.8  1024.7   691.3
## 1995   339.5  1251.9   785.6   440.3
## 1996   641.6  1831.8  1222.4   171.1
## 1997     7.2   229.8    36.3    14.4
## 1998   172.3   754.1    54.8   567.2
## 1999   370.7    23.5   513.0   293.2
## 2000   154.7   725.0  1539.0   228.4
## 2001   274.9   897.7   641.9    48.0
## 2002  1145.5   943.5  2441.1  1141.8
## 2003  1908.7   258.8  4636.6  1094.0
## 2004  5499.4   145.8   752.9  4547.1
## 2005    32.2  2769.0  2241.9  1869.5
## 2006  2056.0  1111.6  2130.7  6534.8
## 2007  1068.8   844.7  1074.5  3477.9
## 2008  2008.2  2838.1  5119.3   194.2
## 2009  1032.9   534.9  3950.7 14537.4
## 2010  3431.1    60.2  4278.2 11344.1
## 2011   747.5  1841.1  4773.5  8888.8
## 2012  6213.9   131.3  4684.3  5046.7
## 2013   682.4   501.7  1720.7  5828.5
## 2014   557.7  1952.7  3132.9  3524.5
## 2015  5732.5  1803.1  1963.2  2751.3
## 2016  2448.9  3626.3  2729.4  6903.8
## 2017  1762.0  3239.2  5112.6  1198.4
## 2018  5057.7  3355.0  1802.6  2702.3
## 2019  2257.7   263.7  3647.2   872.8
## 2020  3472.1 19825.1  5984.0  1034.4
## 2021   473.4   242.5   789.7 13771.2
## 2022  1985.0 14179.1   775.9  2014.8
## 2023

plot(PIB.ts,type="o",col = "black",lwd=2,lty=1)
legend("topright",c("PIB"),lwd=2,lty=1,col=c(1))

plot(pibd,type="o",col = "red",lwd=2,lty=1)
legend("topright",c("Primeira Diferença PIB"),lwd=2,lty=1,col=c(2))

1.1.6 Processo Estocástico

Pode-se definir o processo estocástico como “Aquele que não é determinístico, ou seja, refere-se a uma variável aleatória cujo valor futuro não pode ser previsto com certeza absoluta”. Ou seja, terá um termo de erro $ε$ (uma incerteza). Exemplo: Y = 100A + 0,1B + ε.

Imagine a taxa de câmbio $R\$/US\$$ em cada instante de tempo t entre 17h e 18h neste dia seja aleatório. Pode-se interpretar este fato como uma realização $Z_t(w)$ da variável aleatória $Z_t$, e observar $Z_t(w), 5\lt t\lt 6$. Para se fazer uma previsão às 18h acerca da taxa de câmbio $Z_19(w)$ às 19h, é razoável olhar a evolução total de $Z_t(w)$ entre 17h e 18h. O modelo matemático que descreve esta evolução é chamado de Processo Estocástico.

Podemos definir o Processo Estocástico como “uma coleção de variáveis aleatórias ordenadas no tempo e definidas em um conjunto de pontos T, que pode ser contínuo ou discreto”.

Está-se tentando modelar o caráter aleatório do processo e não uma descrição do tipo causa-efeito como no modelo de regressão. O comportamento pode ser obtido a partir de uma distribuição de probabilidades, e será tanto melhor quanto mais fiel estiver esta distribuição em relação à distribuição verdadeira.

No passeio aleatório, como exemplo, um choque temporário no PIB não se dissipará depois de vários anos e, desta forma, o choque teria efeito de um choque permanente. Se ele se revertesse à tendência, então apenas retirando a tendência o problema estaria resolvido para o uso econométrico.

Da discussão anterior sobre o passeio aleatório, obteve-se $X_t$ como função do valor inicial $X_0=0$, tal que: \[ X_t=X_0+\sum_{j=1}^t{ε_j} \]

Fazendo o valor esperado e a variância de $X_t$ , tem-se a média e a variância dependentes do tempo:

$E(X_t)=\sum_{j=1}^t{E(ε_j)} =tμ$

$Var(X_t)=\sum_{j=1}^t{Var{(ε_j)}} =tσ_ε^2$

O comportamento futuro dos dados no processo estocástico serão descritos pela distribuição de probabilidade conjunta ($F$):

$F(X_1,\ldots,X_t )=P(X_1 \le a_1,\ldots,X_t \le a_t)$

Em que P denota a probabilidade.

Isto posto, pode-se definir a estacionariedade de modo mais explícito, como Estacionariedade Forte. O processo será estacionário forte:

de primeira ordem se: $F(X_{t_1})=F(X_{t_1+k})$
de segunda ordem se: $F(X_{t_1},X_{t_2})=F(X_{t_1+k},X_{t_2+k})$
de ordem n se: $F(X_{t_1},\ldots,X_{t_n})=F(X_{t_1+k},\ldots,X_{t_n+k})$

A Estacionariedade fraca de ordem n pode ser definida quando:

A média e a variância do processo são constantes no tempo; e,
A estrutura de dependência linear depende apenas da distância entre os períodos e diminui com esta distância. Portanto, tem-se Estacionariedade fraca para:

\[ E(X_t) = μ = constante \]

\[ Var(X_t) = \sigma^2 = constante \]

\[ Corr(X_t,X_{t-k}) = \rho(k) \]

2 Função de Autocorrelação e Função de Autocorrelação Parcial

A série temporal dos retornos de um ativo denotada por $\left\{r_{t}\right\}_{t=1}^{T}$ é uma coleção de variáveis aleatórias coletadas ao longo do tempo (definição de processo estocástico). A modelagem econométrica de séries temporais univariadas tem como objetivo capturar a relação linear entre $r_{t}$ e informações disponíveis antes de $t$.

Desta forma, os valores históricos de $r_{t}$ podem ser úteis para modelar o seu comportamento ao longo do tempo. Neste caso, a correlação entre os retornos tem um papel importante. Chamamos tais correlações de autocorrelação e esta é uma ferramenta básica para estudar uma série temporal estacionária.

Neste sentido, este material tem como objetivo contribuir para o entendimento sobre correlação, autocorrelação e autocorrelação parcial.

2.1 Correlação

O coeficiente de correlação de Pearson que mede a dependência linear entre duas variáveis é definido como:

\[ {\rho}_{x,y}=\frac{COV\left(x,y \right)}{\sqrt{Var\left(x \right),Var\left(y \right)}} =\frac{E\left[\left(x-{\mu}_{x} \right) \left(y-{\mu}_{y} \right) \right]}{\sqrt{E{\left(x-{\mu}_{x} \right)}^{2}E{\left(y-{\mu}_{y} \right)}^{2}}} \]

onde ${\mu}_{x}$ e ${\mu}_{y}$ são a média de $x$ e $y$, respectivamente. Algumas propriedades:

$-1\leq {\rho}_{x,y} \leq 1$
${\rho}_{x,y}={\rho}_{y,x}$
${\rho}_{x,y}=0$ quando as duas variáveis não são correlacionada

Existem outras alternativas para o cálculo de correlação proposto por Pearson. Dentre eles temos o $\rho$ de Spearman e o $\tau$ de Kendall.

Quando uma amostra $\left\{(x_{t},y_{t})|t=1,...,T\right\}$ é coletada, a correlação pode ser obtida por seu estimador amostral:

\[ {\hat{\rho}}_{x,y}=\frac{\sum_{t=1}^{T}{ \left({x}_{t}-\overline{x} \right) \left({y}_{t}-\overline{y} \right)}}{ \sqrt{\sum_{t=1}^{T}{{\left({x}_{t}-\overline{x} \right)}^{2}\sum _{t=1}^{T}{{\left({y}_{t}-\overline{y} \right) }^{2}}}}} \]

onde $\overline{x}={\sum_{t=1}^{T}{{x}_{t}}}/{T}$ e $\overline{y}={\sum_{t=1}^{T}{{y}_{t}}}/{T}$ correspondem à média amostral de $x$ e $y$, respectivamente. No gráfico abaixo, mostramos exemplos de variáveis correlacionadas (positiva e negativa).

data(mtcars)

x = sample(1:20,20)+rnorm(10,sd=2)
y = x+rnorm(10,sd=3)
z = (sample(1:20,20)/2)+rnorm(20,sd=5)
df = data.frame(x,y,z)

c1 = hchart(df, "scatter", hcaes(x = x, y = y), color = "black") %>%  
      hc_exporting(enabled = TRUE) 

c2 = hchart(mtcars, "scatter", hcaes(x = wt, y = mpg), color = "red") %>%
      hc_exporting(enabled = TRUE) 
lst = list(c1,c2)

hw_grid(lst, rowheight = 300)  %>% browsable()

2.2 Função de Autocorrelação

A função de autocorrelação é o gráfico da autocorrelação contra a defasagem. Considere uma série temporal de retornos de uma ação $\left\{r_{t}\right\}_{t=1}^{T}$. O coeficiente de correlação entre $r_{t}$ e $r_{t-k}$ é chamado de autocorrelação de k-ésima ordem e é denotadado por:

\[ {\rho}_{k}=\frac {Cov\left({r}_{t},{r}_{t-k} \right)}{\sqrt{Var\left({r}_{t},{r}_{t-k} \right)}} =\frac{Cov\left({r}_{t},{r}_{t-k} \right)}{Var\left({r}_{t} \right)} =\frac{{\gamma}_{k}}{{\gamma}_{0}} \]

onde $Var\left({r}_{t-k}\right)=Var\left({r}_{t}\right)$ porque $r_{t}$ é fracamente estacionário. Além disso, temos:

${\rho}_{0}=1$;
${\rho}_{l}={\rho}_{-l}$;
$-1\leq {\rho}_{l} \leq 1$

Um conjunto de autocorrelações, $\left\{\rho_{k}\right\}$, é chamado de função de autocorrelação de $r_{t}$. Para uma dada amostra de retornos de uma ação, $\left\{r_{t}\right\}_{t=1}^{T}$, suponha que $\overline{r}$ é a média amostral. Então, a autocorrelação amostral de primeira ordem de $r_{t}$ é:

\[ {\hat{\rho}}_{1}=\frac{\sum _{t=2}^{T}{\left({r}_{t}-\overline{r}\right) \left({r}_{t-1}-\overline{r}\right)}}{\sum_{t=1}^{T}{{\left({r}_{t}-\overline{r}\right)}^{2}}} \]

que é um estimador consistente de ${\rho}_{1}$. Em geral, a autocorrelação amostral de k-ésima ordem de $r_{t}$ pode ser definida como:

\[ {\hat{\rho}}_{k}=\frac{\sum_{t=k+1}^{T}{\left({r}_{t}-\overline{r}\right) \left({r}_{t-k}-\overline{r} \right)}}{\sum_{t=1}^{T}{{\left({r}_{t}-\overline{r}\right)}^{2}}} \] para $0\leq k \leq T-1$.

Por exemplo, suponha que você está avaliando uma série temporal qualquer e quer visualizar como as defasagens da série podem impactar seu valor atual (ou seja, se $r_{t}$ é relacionado com $r_{t-k}$ para $k\ge1$). A função de autocorrelação pode ser usada para obter tal informação.

Num primeiro momento, visualize os dados da série para 10 defasagens. Observe que as defasagens se tornam novas colunas e na medida que elas aumentam, incrementa-se as linhas sem observações.

from <- as.Date("1974-01-01")
to <- as.Date("1989-12-31")
days <- seq.Date(from=from,to=to,by="days")
timeseries <- as.zoo((arima.sim(list(order = c(3,0,0), ar = c(0.5,0.15,0.05)), n = 10000)))
lags <- Lag(timeseries, k=1:10)
final <- round(zoo::cbind.zoo(timeseries,lags), 4)
colnames(final) <- c("atual", "lag1", "lag2", "lag3", "lag4", "lag5",
                    "lag6", "lag7", "lag8", "lag9", "lag10")
DT::datatable(final[1:10,], rownames = FALSE, options = list(columnDefs = list(list(className = 'dt-center', targets="_all")), dom = 't'))

Agora, observe a matriz de correlações entre a série temporal e suas defasagens (aqui, apenas 10). O comportamento da correlação é evidenciado pela mudança de cor de azul para azul mais claro, sendo que a cor azul representa correlação positiva e quando vai clareando a correlação diminui. Observamos que há correlação positiva para as primeiras defasagens da série temporal.

hchart(cor(as.data.frame(final), use = "pairwise.complete.obs")) %>% hc_exporting(enabled = TRUE)

Apesar da simples correlação entre os dados nos ajudar a identificar defasagens que poderíam contribuir para o comportamento da série em $t$, precisamos fazer uso de testes estatísticos que verifiquem a significância da relação entre o valor atual e suas defasagens. Neste sentido, a função de autocorrelação tem grande importância.

Abaixo, um exemplo de função de autocorrelação. Observe que há duas linhas horizontais que representam os limites do teste de significância sendo que valores acima ou abaixo da linha são estatisticamente significantes. Neste documento, apresentaremos o teste que é realizado.

# Calcular a autocorrelação 
acf <- stats::acf(timeseries, na.action = na.pass, plot = FALSE, lag.max = 15)

# Gráfico da função de autocorrelação. 
plot(acf, main = "", ylab = "", xlab = "Defasagem")
title("Função de Autocorrelação (FAC)", adj = 0.5, line = 1)

2.3 Função de Autocorrelação Parcial

Em um modelo AR(1), existe uma correlação implícita entre $y_{t}$ e $y_{t-2}$. Isso está presente na FAC, por meio do decaimento exponencial.

Suponha que você quer modelar o retorno no instante $t$, $r_{t}$, como função do retorno imediatamente anterior, $R_{t-1}$. Porém, pode existir também uma correlação implícita entre os retornos em $t$ e $t-2$, mas seu interesse continua sendo verificar se apenas o retorno em $t-1$ é importante para o comportamento do retorno em $t$.

Uma alternativa para filtrar correlações e manter-se apenas a correlação pura entre duas observações é fazer uso da correlação parcial. Formalmente, fazemos uso das seguintes regressões para cada defasagem $j$ de interesse:

\[ {r}_{t}={\phi}_{j,1}{r}_{t-1}+{\phi}_{j,2}{r}_{t-2}+\cdot \cdot \cdot +{\phi}_{j,j}{r}_{t-j}+{\epsilon}_{t} \]

onde ${\epsilon}_{t}$ é um erro e $j=1,2,...$. Em outras palavras, o procedimento faz:

Para $j=1$: regredir $r_{t}$ contra $r_{t-1}$ e obter $\hat{{\phi}}_{1,1}$
Para $j=2$: regredir $r_{t}$ contra $r_{t-1}$ e $r_{t-2}$ e obter $\hat{{\phi}}_{2,1}$ e $\hat{{\phi}}_{2,2}$, mas só interessa $\hat{{\phi}}_{2,2}$
Assim por diante

No caso de uma série temporal, geramos uma função de autocorrelação parcial que será o gráfico da autocorrelação parcial contra possíveis defasagens da própria série temporal, ou seja, o gráfico de $\hat{{\phi}}_{j,j}$ para cada $j$ possível.

Abaixo, gráfico da FACP para a mesma série temporal que usamos para exemplificar a matriz de correlações e a FAC. Observe a diferença entre a função de autocorrelação e função de autocorrelação parcial. Enquanto a primeira mostra que há autocorrelação significante até a defasagem 12, o gráfico da função de autocorrelação parcial define que apenas 3 defasagens da série temporal realmente são importantes para modelar seu valor esperado em $t$.

# Calcular a autocorrelação parcial
pacf <- pacf(timeseries, na.action = na.pass, plot = FALSE, lag.max = 15)

# Gráfico da função de autocorrelação parcial. 
plot(pacf, main = "", ylab = "", xlab = "Defasagem")
title("Função de Autocorrelação Parcial (FACP)", adj = 0.5, line = 1)

Essa avaliação das duas funções é de suma importância para a modelagem de séries temporais lineares dado que contribui para identificar a especificação correta de qual formulação econométrica usar. Neste caso específico, um modelo AR(3), ${r}_{t}={\phi}_{0}+{\phi}_{1}{r}_{t-1}+{\phi}_{2}{r}_{t-2}+{\phi}_{3}{r}_{t-3}+{a}_{t}$, deve ser uma boa especifição. Veremos maiores detalhes sobre esses modelos em próximas aulas.

2.4 Teste de siginificância estatísitica autocorrelação

Para um dado $k$, os resultados da Função de Autocorrelação podem ser testados usando um teste que verifique se:

\[ \begin{aligned} && H_{0}: \rho_{k}=0 \\ && H_{1}: \rho_{k}\neq 0 \end{aligned} \] A estatística do teste será:

\[ t-ratio=\hat{\rho}_{k}\sqrt{T} \]

A decisão do teste será rejeitar $H_{0}$ se $\left| t-ratio \right|>{{Z}_{\alpha}}/{{2}}$, onde ${Z}_{{\alpha}/{2}}$ é o $100(1-{\alpha}/{2})$ percentil de uma distribuição normal.

Exemplo: tenho uma série temporal com $516$ observações e quero testar a signifiância de $\hat{\rho}_{12}=0.13$. Para tanto, basta fazer: $t-ratio=\sqrt{516}*0.13=2.96$ que é maior em módulo do que o valor crítico de $5\%$ que é $1.96$. Conclusão: rejeitamos $H_{0}$.

2.5 Testar conjutamente várias autocorrelações

As estatísticas $\hat{\rho}_{1}$, $\hat{\rho}_{2}$, … são chamadas de FAC amostral de uma série temporal. Elas têm um papel importante na análise de séries temporais lineares. Em muitas aplicações estamos interessados em testar se várias defasagens da função de autocorrelação são iguais a 0. Box - 1970 propos um teste, conhecido como Ljung-Box, que foi modificado por Ljung - 1978 para melhorar o poder do teste em amostras finitas e é definido como:

\[ Q=T(T+2)\sum_{j=1}^{n}{\frac{{\hat{\rho}}_{j}^{2}}{T-j}} \] que segue uma distribuição ${\chi}_{n}^{2}$ que indica $n$ graus de liberdade. Sejam as hipóteses nulda, dada por $H_{0}:\sum_{j=1}^{n}{\rho_{j}=0}$, e alternativa, dada por $H_{0}:\sum_{j=1}^{n}{\rho_{j}\neq0}$, então podemos testar conjuntamente se $j$ defasagens da função de autocorrelação são iguais a $0$.

2.6 FAC e FACP na prática

Mas antes de chegar nas Funções de autocorrelação e função de autocorrelação parcial quero deixair claro com exemplos o que abordaremos a seguir:

knitr::include_graphics('./decaimento.png')

Este é um exemplo de grafico com decaimento ate zero de forma geometrica! ☝️

knitr::include_graphics('./senoide.png')

Este é um exemplo de senóide amortecida! ☝️

knitr::include_graphics('./facp.png')

E este é um exemplo de quando falarmos que havera um decaimento brusco com algum (s) spike/prego significante (s)!☝️

2.6.1 Modelo AR

dataframe <- read_excel(path = "AR.xlsx", col_names=TRUE)

AR1 <- ts(dataframe, start = 2010, frequency = 12)
AR1

##         Jan    Feb    Mar    Apr    May    Jun    Jul    Aug    Sep
## 2010  0.656  1.057 -1.750 -0.489 -2.861 -2.227 -2.014 -3.773 -3.333
## 2011 -1.801 -0.538 -0.292 -0.444  1.648  2.183 -0.253 -1.069 -2.092
## 2012  3.098  4.853  4.649  4.821  4.441  5.496  3.892  4.290  3.746
## 2013  2.144  1.252 -0.006  0.412  1.958  1.883  0.344 -0.708 -1.852
## 2014 -1.799 -1.698                                                 
##         Oct    Nov    Dec
## 2010 -0.626 -0.731 -0.549
## 2011 -1.993 -1.187  1.394
## 2012  3.723  1.111  3.480
## 2013 -2.318 -2.300 -0.937
## 2014

ggtsdisplay(AR1)

Pela analise da FAC, nós percebemos que ela tem um padrao de decaimento até o zero de forma geometrica. Já a FACP apresenta um decaimento brusco e somente com 1 spike (prego) significante assim nos indicado que devemos trabalhar com um modelo AR(1)

dataframe <- read_excel(path = "AR_.xlsx", col_names=TRUE)

AR2 <- ts(dataframe, start = 2010, frequency = 12)
AR2

##       Jan  Feb  Mar  Apr  May  Jun  Jul  Aug  Sep  Oct  Nov  Dec
## 2010 5.16 5.34 4.99 4.70 4.90 4.92 5.11 4.99 5.19 5.05 4.84 5.08
## 2011 4.86 4.95 4.73 4.90 4.94 4.83 4.99 5.19 5.34 5.13 4.90 4.33
## 2012 4.70 5.42 5.69 5.61 4.83 4.82 4.97 4.99 4.81 5.13 5.28 5.07
## 2013 4.88 4.81 4.48 4.62 4.94 5.48 5.09 4.55 4.56 5.19 5.38 5.03
## 2014 4.91 4.55 4.47 4.98 5.27 5.08 4.79 4.71 4.77 5.31 5.37 4.86
## 2015 4.73 4.73 5.20 5.00 5.04 5.10 4.78 4.82 4.50 4.70 4.74 4.78
## 2016 4.86 5.02 5.20 5.00 4.80 5.32 5.20 5.21

ggtsdisplay(AR2)

Pela analise da FAC, nós percebemos que ela tem um padrao de senóide amortecida (outra forma de se apresentar, mas que ao fim é igual ao exemplo anterior onde se apresentava um padrao de decaimento até o zero de forma geometrica, a unica diferença é que aqui os valores alternam entre positivos e negativos, assim resultanto nessas diversas curvas). Já a FACP apresenta um decaimento brusco, mas agora com 2 spikes (pregos) significanetes, assim nos indicando que devemos trabalhar com um modelo AR(2)

2.6.2 Modelo MA

dataframe <- read_excel(path = "MA.xlsx", col_names=TRUE)

MA <- ts(dataframe, start = 2010, frequency = 12)
MA

##         Jan    Feb    Mar    Apr    May    Jun    Jul    Aug    Sep
## 2010  0.811 -0.028  0.159  0.663 -2.299 -0.494  1.498 -1.408  0.032
## 2011 -0.316 -1.513  0.545  0.071  0.447 -0.287  0.570 -0.336 -2.074
## 2012 -1.902  2.908 -3.726  1.376  2.261 -2.789 -0.220  1.038  0.730
## 2013 -0.974  0.524  1.093 -0.118 -2.390  1.765 -0.038  0.249 -0.460
## 2014 -0.497  0.288                                                 
##         Oct    Nov    Dec
## 2010  0.487  0.888  0.458
## 2011  0.748  0.858  0.248
## 2012  0.664  1.166 -1.548
## 2013  0.894 -0.775 -0.031
## 2014

ggtsdisplay(MA)

Agora o comportamento da FAC e FACP estão invertidos! Ou seja estamos trabalhando com Médias Moveis. Portanto sua FAC apresenta padrão de decaimento brusco ate zero, nesse caso temos 1 spike (prego) significante, portanto trabalhariamos com um modelo MA(1). E a FACP de um modelo MA apresenta padrao de decaimento até zero de maneira geometria ou formato senóide amortecido. Nesse caso especifico temos a primeira possibilidade mesmo (Mas de cabeça pra baixo 🙃)

2.6.3 Modelo ARMA

dataframe <- read_excel(path = "ARMA.xlsx", col_names=TRUE)

ARMA <- ts(dataframe, start = 2000, frequency = 12)
ARMA

##         Jan    Feb    Mar    Apr    May    Jun    Jul    Aug    Sep
## 2000  0.466  0.069  1.158  0.187  0.483  0.920  0.754  0.218  0.810
## 2001  0.705  1.009  0.147  1.336 -0.793 -0.484 -1.932 -1.333 -0.752
## 2002  1.336  0.904  1.836 -0.259  1.876  3.688  1.988  2.241  2.993
## 2003  1.929  0.874  1.556  0.267  1.496 -0.207  2.265  2.819  3.425
## 2004 -0.047  1.100                                                 
##         Oct    Nov    Dec
## 2000  1.250  0.949  0.913
## 2001  0.092 -0.562  1.521
## 2002  3.515  3.185  4.370
## 2003  1.072  1.150  1.340
## 2004

ggtsdisplay(ARMA)

Agora para os Modelos ARMA, tanto a FAC quanto FACP terão um padrão de decaimento até o zero de maneira exponencial ou formato senóide amortecido. Em nosso caso, os dois modelos apresentam um padrão de senóide amortecida. E como vamos saber com quantas defasagens vamos trabalhar? É simples, A FAC nos da as informações de MA e a FACP as informações do AR. Devemos olhar quantos spikes REALMENTE diferenciados temos e assim decidir o numero de defasagens. Em nosso caso temos um Modelo ARMA(1,1).

Aula - FAC, FACP e a Estacionariedade

Daniel, Lucas e Luiz

18 junho 2023